##########################################################################################
# PhyOP pipeline Makefile



#submit job to cluster 
CMD_SUBMIT?=qrsh -now n -cwd -q short_jobs.q -p -5
#CMD_SUBMIT=qrsh -now n -cwd -q bc1.q,bc2.q,highmem.q -p -5
#CMD_SUBMIT=./fake_queue.pl
GENOME_ALIGN_MAK=genome_align.mak

SEQ_BATCH_SIZE?=200



# set per project
#PROJ_NAME=homo_monodelphis
#SPECIES1=homo
#SPECIES2=monodelphis
#REGEX1=(ENS[^ \t]+)
#REGEX2=(mono[^ \t]+)
#DIR_PIPELINE_TMP=/net/cpp-data/scratch/Leo/genome_temp/monodelphis_project/pipeline_ds_temp/
#DIR_PROJECT=/net/cpp-group/Leo/monodelphis

# shouldn't need to be changed often
DIR_BIN?=/net/cpp-group/Leo/bin



#
#    Calculated
#
DIR_SEQ?=$(DIR_PROJECT)/sequences
DIR_PROJECTBIN?=$(DIR_PROJECT)/bin
SPECIES12=$(SPECIES1)_$(SPECIES2)
DIR_ALIGNRES?=$(DIR_PROJECT)/align_results

# for convenience
DIR_SEQ1=$(DIR_SEQ)/$(SPECIES1)
DIR_SEQ2=$(DIR_SEQ)/$(SPECIES2)
DIR_ALIGNED1=$(DIR_ALIGNRES)/$(SPECIES1)
DIR_ALIGNED2=$(DIR_ALIGNRES)/$(SPECIES2)
DIR_ALIGNED12=$(DIR_ALIGNRES)/$(SPECIES12)


#
#    Create directories if necessary
#
DIR_OUTPUT?=$(DIR_PROJECT)/orthologs/$(SPECIES12)
DIR_KAKSRES?=$(DIR_PROJECT)/kaks_results/$(SPECIES12)

# pipeline
DIR_P_SUMMARY=$(DIR_PIPELINE_TMP)/summary
DIR_P_ERRORS=$(DIR_PIPELINE_TMP)/errors
DIR_P_KAKS=$(DIR_PIPELINE_TMP)/kaks
DIR_P_KAKS_RESULTS=$(DIR_PIPELINE_TMP)/kaks_results
DIR_P_TREE=$(DIR_PIPELINE_TMP)/tree
DIR_P_GCNT=$(DIR_PIPELINE_TMP)/gene_counts
DIR_P_FILTERED_KAKS=$(DIR_PIPELINE_TMP)/filtered_kaks
DIR_P_BLASTRES=$(DIR_PIPELINE_TMP)/blast_results
DIR_P_SEQ1=$(DIR_PIPELINE_TMP)/sequences/$(SPECIES1)
DIR_P_SEQ2=$(DIR_PIPELINE_TMP)/sequences/$(SPECIES2)
DIR_P_SEQ12=$(DIR_PIPELINE_TMP)/sequences/$(SPECIES12)

# for convenience
DIR_ALIGNED1=$(DIR_ALIGNRES)/$(SPECIES1)
DIR_ALIGNED2=$(DIR_ALIGNRES)/$(SPECIES2)
DIR_ALIGNED12=$(DIR_ALIGNRES)/$(SPECIES12)

# output
DIR_OUTPUT_GENE_LISTS=$(DIR_OUTPUT)/gene_lists
DIR_OUTPUT_INTERPRET_TREE=$(DIR_OUTPUT)/interpret_tree



#
# Make in three parts. Second part depends on first etc
#

all: variable_prerequisites do_fifth_part




.DELETE_ON_ERROR:
.PHONY: all are_variables_defined first_part print_variables
.SECONDARY :	


##########################################################################################
#
#   Functions


#
#	prints only last directory+file
#
#	$(call dir_file,FILE_NAME)

dir_file=$(notdir $(patsubst %/,%,$(dir $(1))))/$(notdir $(1))

#
#	create file unless already exists
#
#	$(call create_unless_exists,FILE_NAME)

create_unless_exists= @if [[ ! -f $(1) ]]; then \
echo "	$(call dir_file,$(1)) is not present. Create empty placeholder"; touch $(1); fi


#
#	blastall if both database and fasta non-zero
#
#	$(call blastall_if_non_zero,NAME,DATABASE,FA_FILE,RESULT)
blastall_if_non_zero=@if [[ -s $(2) && -s $(3) ]]; then \
	echo "	Blastall $(1) submitted"; \
	@$(CMD_SUBMIT) -N "blst_$(1)" "$(DIR_BLASTBIN)/blastall -M BLOSUM80 -p blastp -d $(2) -a 2 -i $(3) -e 1e-5 -F F -b 100000 -v 100000 -a 2 -X 150 -m 7 -z $(EFFECTIVE_LEN_BLAST_DB) | $(DIR_BIN)/parse_blastall_output.pl --reciprocal > $(4)" ;\
	echo "	Blastall $(1) finished"; \
else \
	rm -f $(4); \
	touch $(4); \
	echo "		Blastall $(1) skipped" ; \
fi

#
#	consolidate_non_empty_files
#
#	$(call consolidate_non_empty_files,DIR,WILDCARD,RESULT_FILE)
consolidate_non_empty_files= \
        @echo "	Consolidate to $(3)"; \
	rm -f $(3); \
        find $(1) -name "$(2)" -empty -print0 \
	| xargs -0 -i rm {} ; \
        find $(1) -name "$(2)"  -print0 \
        | xargs -0 -i bash -c "echo '{}' >> $(3); cat '{}' >> $(3); echo >> $(3)"; \
        find $(1) -name "$(2)"  -print0 | xargs -0 -i rm {}
#
#
#
#
##########################################################################################



##########################################################################################
#
#   Make directories
#
ALL_DIR= $(DIR_OUTPUT) \
         $(DIR_KAKSRES) \
         $(DIR_P_SUMMARY) \
         $(DIR_P_ERRORS) \
         $(DIR_P_KAKS) \
         $(DIR_P_TREE) \
         $(DIR_P_GCNT) \
         $(DIR_P_FILTERED_KAKS) \
         $(DIR_P_BLASTRES) \
         $(DIR_P_SEQ1) \
         $(DIR_P_SEQ2) \
         $(DIR_P_SEQ12) \
         $(DIR_P_KAKS_RESULTS) \
         $(DIR_OUTPUT_INTERPRET_TREE) \
         $(DIR_OUTPUT_GENE_LISTS)


makedirectories: $(ALL_DIR)

$(ALL_DIR):
	@echo "	Make $@ directory"
	@mkdir -p $@
#
#
##########################################################################################




##########################################################################################
#
#	get cdna sequences for species pair
#
#	
get_cdna_sequences: $(DIR_P_SEQ12)/cdna.fa
	@echo "	Get cDNA sequences"

$(DIR_P_SEQ1)/cdna.fa: $(DIR_SEQ1)/cdna.mapped.fa
	@cp $< $@
$(DIR_P_SEQ2)/cdna.fa: $(DIR_SEQ2)/cdna.mapped.fa
	@cp $< $@

$(DIR_P_SEQ12)/cdna.fa: $(DIR_SEQ1)/cdna.mapped.fa $(DIR_SEQ2)/cdna.mapped.fa
	@cat $^ > $@

#
#
##########################################################################################

##########################################################################################
#
#	get identifier maps for species pair
#
#	

get_id_maps: $(DIR_P_SEQ12)/cdna_hash_to_prot_id.map \
             $(DIR_P_SEQ12)/cdna_hash_to_gene_id_species.map \
	     $(DIR_OUTPUT)/prot_to_gene_ids.map
	@echo "	Get Identifier maps"

$(DIR_P_SEQ12)/cdna_hash_to_gene_id_species.map:  $(DIR_SEQ1)/cdna_hash_to_gene_id_species.map \
                                                  $(DIR_SEQ2)/cdna_hash_to_gene_id_species.map
	@cat $^ > $@

$(DIR_P_SEQ12)/cdna_hash_to_prot_id.map:  $(DIR_SEQ1)/cdna_hash_to_prot_id.map \
                                          $(DIR_SEQ2)/cdna_hash_to_prot_id.map
	@cat $^ > $@


$(DIR_OUTPUT)/prot_to_gene_ids.map:  $(DIR_SEQ1)/prot_to_gene_ids.map \
                                      $(DIR_SEQ2)/prot_to_gene_ids.map
	@cat $^ > $@





#
#
##########################################################################################


##########################################################################################
#
#   Check variables defined
#
#
are_variables_defined:
ifeq (undefined, $(origin DIR_PIPELINE_TMP))
	@echo "	DIR_PIPELINE_TMP is not undefined"; exit 255
endif
ifeq (undefined, $(origin SPECIES1))
	@echo "	SPECIES1 is not undefined"; exit 255
endif
ifeq (undefined, $(origin SPECIES2))
	@echo "	SPECIES2 is not undefined"; exit 255
endif
ifeq (undefined, $(origin REGEX1))
	@echo "	REGEX1 is not undefined"; exit 255
endif
ifeq (undefined, $(origin REGEX2))
	@echo "	REGEX2 is not undefined"; exit 255
endif
ifeq (undefined, $(origin DIR_PROJECT))
	@echo "	DIR_PROJECT is not undefined"; exit 255
endif
ifeq (undefined, $(origin PROJ_NAME))
	@echo "	PROJ_NAME is not undefined"; exit 255
endif
#
##########################################################################################







##########################################################################################
#
#   Print variables
#
#

print_variables:
	@echo
	@echo "	****************************************************************************"
	@echo
	@echo "Variables:"
	@echo "	PROJ_NAME	=$(PROJ_NAME)         "
	@echo "	SPECIES1	=$(SPECIES1)          "
	@echo "	SPECIES2	=$(SPECIES2)          "
	@echo "	REGEX1		=$(REGEX1)            "
	@echo "	REGEX2		=$(REGEX2)            "
	@echo "	DIR_PIPELINE_TMP=$(DIR_PIPELINE_TMP)  "
	@echo "	DIR_PROJECT	=$(DIR_PROJECT)          "
	@echo "	DIR_BIN		=$(DIR_BIN)           "
	@echo "	DIR_PROJECTBIN	=$(DIR_PROJECTBIN)    "
	@echo "	DIR_SEQ		=$(DIR_SEQ)           "
	@echo "	DIR_ALIGNRES	=$(DIR_ALIGNRES)      "
	@echo "	DIR_KAKSRES     =$(DIR_KAKSRES)       "
	@echo "	SPECIES12	=$(SPECIES12)      "
	@echo "	DIR_P_SUMMARY	=$(DIR_P_SUMMARY)     "
	@echo "	DIR_P_ERRORS	=$(DIR_P_ERRORS)      "
	@echo "	DIR_P_SEQ1	=$(DIR_P_SEQ1)        "
	@echo "	DIR_P_SEQ2	=$(DIR_P_SEQ2)        "
	@echo "	DIR_P_NEWNEW1	=$(DIR_P_NEWNEW1)     "
	@echo "	DIR_P_NEWNEW2	=$(DIR_P_NEWNEW2)     "
	@echo "	DIR_P_OLDOLD1	=$(DIR_P_OLDOLD1)     "
	@echo "	DIR_P_OLDOLD2	=$(DIR_P_OLDOLD2)     "
	@echo "	DIR_P_OLDNEW1	=$(DIR_P_OLDNEW1)     "
	@echo "	DIR_P_OLDNEW2	=$(DIR_P_OLDNEW2)     "
	@echo "	DIR_P_NEWNEW12	=$(DIR_P_NEWNEW12)    "
	@echo "	DIR_P_NEW2OLD1	=$(DIR_P_NEW2OLD1)    "
	@echo "	DIR_P_NEW1OLD2	=$(DIR_P_NEW1OLD2)    "
	@echo "	DIR_P_OLD1OLD2	=$(DIR_P_OLD1OLD2)    "
	@echo "	DIR_SEQ1	=$(DIR_SEQ1)          "
	@echo "	DIR_SEQ2	=$(DIR_SEQ2)          "
	@echo "	DIR_ALIGNED1	=$(DIR_ALIGNED1)      "
	@echo "	DIR_ALIGNED2	=$(DIR_ALIGNED2)      "
	@echo "	DIR_ALIGNED12	=$(DIR_ALIGNED12)     "
	@echo " DIR_P_KAKS      =$(DIR_P_KAKS)        "
	@echo " DIR_P_TREE      =$(DIR_P_TRE)         "
	@echo " DIR_P_GCNT      =$(DIR_P_GCNT)        "
	@echo " DIR_P_FILTERED_KAKS=$(DIR_P_FILTERED_KAKS)"
	@echo
	@echo "	****************************************************************************"
	@echo

#
#
##########################################################################################


.DELETE_ON_ERROR:
.PHONY: variable_prerequisites

MYMAKE_FILE= -f $(CURDIR)/makefile.mak


#
# make sure all variable defined
#
variable_prerequisites:
ifeq (undefined, $(origin DIR_PIPELINE_TMP))
	@echo "	DIR_PIPELINE_TMP is not undefined"; exit 255
endif
ifeq (undefined, $(origin SPECIES1))
	@echo "	SPECIES1 is not undefined"; exit 255
endif
ifeq (undefined, $(origin SPECIES2))
	@echo "	SPECIES2 is not undefined"; exit 255
endif
ifeq (undefined, $(origin REGEX1))
	@echo "	REGEX1 is not undefined"; exit 255
endif
ifeq (undefined, $(origin REGEX2))
	@echo "	REGEX2 is not undefined"; exit 255
endif
ifeq (undefined, $(origin DIR_PROJECT))
	@echo "	DIR_PROJECT is not undefined"; exit 255
endif
ifeq (undefined, $(origin PROJ_NAME))
	@echo "	PROJ_NAME is not undefined"; exit 255
endif













##########################################################################################
#
# Cleanup
#
clean:
#
#
#
##########################################################################################









##########################################################################################
#
#   First Part:
#
#           Get all pairs
#
#
do_first_part: makedirectories get_id_maps
#_________________________________________________________________________________________
#
# Run Genome genome align if shared results older than sequences (or identifiers)
#	assume that shared results is always the oldest of species1/species2/shared results
#	because species1/species2 could have been updated by some other species pair analysis
#
$(DIR_ALIGNED12)/curr.blastall_results: $(DIR_SEQ1)/cdna.ids  $(DIR_SEQ2)/cdna.ids 
	@echo "	Genome align necessary..."
	@echo "	****************************************************************************"
	@echo
	@$(MAKE) -f $(GENOME_ALIGN_MAK) all
	@echo "- Genome align Finished"



##########################################################################################
#
#   Second Part:
#
#		Prepare files for Kaks
#
#
do_second_part: do_first_part
	@echo "- Second Part"
	@echo "	Prepare files for Kaks"
	@echo "	Combine alignment results and the then split up for kaks"
	@echo "	****************************************************************************"
	@echo
	@$(MAKE) $(MYMAKE_FILE)  second_part
	@echo "- Second Part Finished"
	@echo

.SECONDARY : $(DIR_P_BLASTRES)/all.align.aln


	

#_________________________________________________________________________________________
#
# Combine all results
#
$(DIR_ALIGNED1)/curr.blastall_results: $(DIR_ALIGNED12)/curr.blastall_results

$(DIR_ALIGNED2)/curr.blastall_results: $(DIR_ALIGNED1)/curr.blastall_results

$(DIR_P_BLASTRES)/curr.blastall_results: $(DIR_ALIGNED2)/curr.blastall_results
	@echo "	Combine all alignment results"
	@cat $(DIR_ALIGNED12)/curr.blastall_results \
            $(DIR_ALIGNED1)/curr.blastall_results  \
            $(DIR_ALIGNED2)/curr.blastall_results  > $@
	@echo "	All alignment results combined"
	@echo



#_________________________________________________________________________________________
#
# Find all alignments without previous kaks results
#
$(DIR_P_BLASTRES)/align.aln.without_kaks: $(DIR_P_BLASTRES)/curr.blastall_results
	@echo "	Previous kaks results can be reused"
	$(call create_unless_exists,$(DIR_KAKSRES)/curr.kaks_results)
	@$(DIR_BIN)/seq_pairs_matching \
		--file1in $(DIR_KAKSRES)/curr.kaks_results \
		--file2in $< \
		--shared_in_full $(DIR_P_KAKS_RESULTS)/matching.previous_kaks_res \
		--file2_only_in_full $@ \
		--counts $(DIR_P_SUMMARY)/blast_res.previous_results.counts \
		--verbose
	@echo

# this is a secondary product of the above procedure
$(DIR_P_KAKS_RESULTS)/matching.previous_kaks_res: $(DIR_P_BLASTRES)/align.aln.without_kaks

#_________________________________________________________________________________________
#
# break up alignment results into chunks of 1000 in preparation of kaks
#
$(DIR_P_KAKS)/align.0.aln: $(DIR_P_BLASTRES)/align.aln.without_kaks
	@echo "	Break into chunks of 1000..."
	touch $@
	@$(DIR_BIN)/break_into_chunks.pl \
	--every 2000 \
	--base_name "$(DIR_P_KAKS)/align.aln" \
	< $<
	@echo `find $(DIR_P_KAKS) -name align.*.aln | wc -l` aligned files for kaks > $(DIR_P_SUMMARY)/kaks.aligned_files.counts
	@date > $(DIR_P_SUMMARY)/kaks.kaks_start.date
	@echo



second_part: $(DIR_P_KAKS)/align.0.aln
#
# End of second part
#
#-----------------------------------------------------------------------------------------










##########################################################################################
#
# Third Part:
#
#
# Kaks
KAKS_INPUT_FILES=$(wildcard $(DIR_P_KAKS)/align.*.aln)
KAKS_RES_FILES=$(KAKS_INPUT_FILES:%.aln=%.kaks_res)


do_third_part: do_second_part get_cdna_sequences
	@echo "- Third Part"
	@echo "	Kaks calculations, recombine"
	@echo "	****************************************************************************"
	@echo
	@$(MAKE) $(MYMAKE_FILE) third_part
	@echo "- Third Part Finished"
	@echo

.SECONDARY : $(KAKS_RES_FILES) $(DIR_KAKSRES)/curr.kaks_results



##########################################################################################
#
#
#	Kaks if aln file is not empty
#
#	$(call kaks_if_non_zero,NAME,ALN,RESULT)
kaks_if_non_zero=@if [[ -s $(2) ]]; then \
	echo "	Kaks $(1) submitted"; \
	$(CMD_SUBMIT) -N "kaks_$(*F)" "$(DIR_BIN)/seq_pairs_kaks \
                --cdna_sequences $(DIR_P_SEQ12)/cdna.fa \
                --temp_directory $*.tmp \
                --err_log $*.kaks.errors \
                --paml_err_log $*.codeml.errors \
                < $(2) \
                > $(3)" ;\
	echo "	Kaks $(1) finished"; \
else \
	rm -f $(3); \
	touch $(3); \
	echo "		Kaks $(1) skipped" ; \
fi
#
#
#
##########################################################################################


#_________________________________________________________________________________________
#
# Run Kaks
#
%.kaks_res: %.aln
	@echo "	Kaks start"
	$(call kaks_if_non_zero,kaks_$(*F),$<,$@)
	@echo "	Kaks finish for $@"
	@touch $@.kaks_finished



#_________________________________________________________________________________________
#
# Concatenate
#
$(DIR_P_KAKS_RESULTS)/kaks.res.all: $(DIR_P_KAKS_RESULTS)/matching.previous_kaks_res $(KAKS_RES_FILES) 
	@echo "	Concatenate all Kaks results"
	@rm -f $@
	@touch $@
	@find $(DIR_P_KAKS) -name "align.*.kaks_res" -print0 | xargs -0 cat >> $@ 
	@cat $< >> $@
	@echo "	finish"
	@date > $(DIR_P_SUMMARY)/kaks.kaks_finish.date
	@echo


#_________________________________________________________________________________________
#
# Save to prev results
#
$(DIR_KAKSRES)/curr.kaks_results: $(DIR_P_KAKS_RESULTS)/kaks.res.all
	@echo " Save with previous kaks results"
	$(call create_unless_exists,$(DIR_KAKSRES)/curr.kaks_results)
	@# create to temp file first in case fails
	@$(DIR_BIN)/seq_pairs_matching \
		--file1in $@ \
		--file2in $< \
		--merged_in_full $@.tmp \
		--verbose
	@# create prev file if necessary and backup to .prev
	@$(call create_unless_exists,$(DIR_KAKSRES)/kaks_results.prev)
	@find $(DIR_KAKSRES) -name "*.prev" -exec mv '{}' '{}'.prev ';'
	@mv $(DIR_KAKSRES)/curr.kaks_results $(DIR_KAKSRES)/kaks_results.prev
	@# write out result file
	@mv $@.tmp $@
	@date > $(DIR_P_SUMMARY)/kaks.$(SPECIES12).kaks_finish.date
	@echo "	Kaks results cached"



$(DIR_P_KAKS_RESULTS)/prot_id.kaks_results: $(DIR_P_KAKS_RESULTS)/kaks.res.all $(DIR_KAKSRES)/curr.kaks_results
	@echo "	Convert from hashes to protein identifiers"
	@$(DIR_BIN)/seq_pairs_substitute_identifiers \
		--input_file $< \
		--map_identifiers $(DIR_P_SEQ12)/cdna_hash_to_prot_id.map \
		--Unmatched_identifiers $(DIR_P_ERRORS)/kaks.cdna_hash_to_prot_ids.missing_ids.errors \
		--err_log $(DIR_P_ERRORS)/kaks.cdna_hash_to_prot_ids.errors \
		--verbose \
		> $@
	@echo "	Finished converting from hashes to protein identifiers"
	


#-----------------------------------------------------------------------------------------
#
#
# Filter by dS using 60% of longer
#
$(DIR_P_FILTERED_KAKS)/ks.res.longer.filtered: $(DIR_P_KAKS_RESULTS)/prot_id.kaks_results
	@echo "	Filter dS results: 60% coverage of longer sequence ..."
	@$(DIR_BIN)/seq_pairs_filter_kaks_results \
		--err_log $(DIR_P_ERRORS)/kaks.filter.longer.errors \
		--coverage_residues 0.6 \
		--coverage_fragment 0.0 \
		--verbose \
		--output_ds_only \
		--results_summary $(DIR_P_SUMMARY)/kaks.filter.longer.counts \
		--input_file $< \
		> $@
	@echo "	Filter Kaks finished"
	@echo


#-----------------------------------------------------------------------------------------
#
#
# Cluster by filtered by longer first
#
$(DIR_P_FILTERED_KAKS)/ks.res.clustered.by.longer.ids: $(DIR_P_FILTERED_KAKS)/ks.res.longer.filtered
	@echo "	First cluster dS using data filtered by longer..."
	@$(DIR_BIN)/seq_pairs_cluster_by_single_linkage \
		--identifiers_only \
		--err_log $(DIR_P_ERRORS)/kaks.cluster.longer.errors \
		--verbose \
		--input_file $< \
		--minimum_cluster_size 10 \
		--results_summary $(DIR_P_SUMMARY)/kaks.cluster.longer.counts \
		> $@
	@echo "	Cluster KaKs results using data filtered by longer finished"
	@echo

#-----------------------------------------------------------------------------------------
#
#
# Filter by dS using 75% of shorter and 65% of longer
#
$(DIR_OUTPUT)/ds_res.both_filtered: $(DIR_P_KAKS_RESULTS)/prot_id.kaks_results $(DIR_P_FILTERED_KAKS)/ks.res.clustered.by.longer.ids
	@echo "	Filter dS : 65% coverage of longer, 75% of shorter sequence ..."
	@$(DIR_BIN)/seq_pairs_filter_kaks_results \
		--err_log $(DIR_P_ERRORS)/kaks.filter.both.errors \
		--coverage_residues 0.6 \
		--coverage_fragment 0.75 \
		--verbose \
		--output_ds_only \
		--results_summary $(DIR_P_SUMMARY)/kaks.filter.both.counts \
		--input_file $< \
		> $@
	@echo "	Filter Kaks finished"
	@echo



#-----------------------------------------------------------------------------------------
#
#
# Save matching kaks
#
$(DIR_OUTPUT)/dnds_res.both_filtered: $(DIR_P_KAKS_RESULTS)/prot_id.kaks_results $(DIR_OUTPUT)/ds_res.both_filtered
	@echo "	dNdS results: 65% coverage of longer, 75% of shorter sequence ..."
	@$(DIR_BIN)/seq_pairs_matching \
		--err_log $(DIR_P_ERRORS)/kaks.filter.matching_both.errors \
		--file1in $< \
		--file2in $(DIR_OUTPUT)/ds_res.both_filtered \
		--shared_in_full $@ \
		--verbose
	@echo "	Done"
	@echo



#-----------------------------------------------------------------------------------------
#
#
# Cluster filtered by both 
#
$(DIR_P_FILTERED_KAKS)/ks.res.clustered: $(DIR_P_FILTERED_KAKS)/ks.res.clustered.by.longer.ids $(DIR_OUTPUT)/ds_res.both_filtered
	@echo "	Parachute into clusters data filtered by shorter ..."
	@$(DIR_BIN)/seq_pairs_cluster_by_single_linkage \
		--predefined_clusters $< \
		--err_log $(DIR_P_ERRORS)/kaks.cluster.both.errors \
		--verbose \
		--results_summary $(DIR_P_SUMMARY)/kaks.cluster.both.counts \
		< $(DIR_OUTPUT)/ds_res.both_filtered \
		> $@
	@echo "	Cluster dS finished"
	@echo


#-----------------------------------------------------------------------------------------
#
#
# Split ds clusters
#
$(DIR_P_TREE)/clusters.0.ks: $(DIR_P_FILTERED_KAKS)/ks.res.clustered
	@echo "	Split dS clusters fragments"
	@$(DIR_BIN)/seq_pairs_split_clusters_by_size.pl \
		--path $(DIR_P_TREE)/clusters.ks \
		< $<
	@echo "	Split dS clusters finish"
	@echo `find $(DIR_P_TREE) -name clusters.*.ks | wc -l` files for kitsch > $(DIR_P_SUMMARY)/tree.ds_clusters.files.counts
	@date > $(DIR_P_SUMMARY)/kaks.upgma_kitsch_start.date
	@echo







third_part: $(DIR_P_TREE)/clusters.0.ks
#
# End of third part
#
#-----------------------------------------------------------------------------------------














##########################################################################################
#
# Fourth Part:
#
#
do_fourth_part: do_third_part
	@echo "- Fourth Part"
	@echo "	Cluster genes with UPGMA and extract large sub-branches for optimisation"
	@echo "	****************************************************************************"
	@echo
	@$(MAKE) $(MYMAKE_FILE)  fourth_part
	@echo "- Fourth Part Finished"
	@echo
	@echo

KS_FILES=$(wildcard $(DIR_P_TREE)/*.ks)
AFTER_FILES=$(KS_FILES:%.ks=%.after)
UPGMA_FILES=$(KS_FILES:%.ks=%.upgma)
.SECONDARY : $(UPGMA_FILES)
# $(UPGMA_FILES)

%.upgma: %.ks
	@echo "	upgma submitted"
	@$(DIR_BIN)/tree_upgma \
		--matrix_output \
		--input_format p \
		--output_format b \
		--header \
		--err_log $(DIR_P_ERRORS)/upgma.$(*F).errors \
		< $*.ks > $*.upgma
	@echo "	upgma finished"
# --verbose

%.after: %.upgma
	@echo "	extract_sub_trees submitted"
	@$(DIR_BIN)/tree_extract_sub_trees_with_matrix \
		--unchanged_tree_file $*.nh.unchanged \
		--err_log $(DIR_P_ERRORS)/extract_sub_tree.$(*F).errors \
		--tree_files_path $*.orig_sub_tree \
		--sub_matrix_files_path $*.sub_matrix \
		--description "$*.after" \
		--input_format b \
		--output_format b < $*.upgma > $*.after
	@echo "	extract_sub_trees finished"
# --verbose



fourth_part: $(AFTER_FILES)
#
# End of fourth part
#
#-----------------------------------------------------------------------------------------



##########################################################################################
#
# Fifth Part:
#
#
do_fifth_part: do_fourth_part
	@echo "- Fifth Part"
	@echo "	Use Kitsch to optimise all large sub-branches"
	@echo "	****************************************************************************"
	@echo
	@$(MAKE) $(MYMAKE_FILE)  fifth_part
	@echo "- Fifth Part Finished"
	@echo
	@echo
SUB_MATRIX_FILES=$(wildcard $(DIR_P_TREE)/*.sub_matrix)
KITSCH_RESULTS_FILES=$(SUB_MATRIX_FILES:%.sub_matrix=%.kitsch)
.SECONDARY : $(KITSCH_RESULTS_FILES) $(DIR_OUTPUT_INTERPRET_TREE)/ds.tree $(DIR_P_TREE)/rejoin.errors

# Create Kitsch from sub matrix
%.kitsch: %.sub_matrix
	@echo "	kitsch submitted"
	@$(CMD_SUBMIT) -N "kitsch_$(*F)" $(DIR_BIN)/tree_kitsch \
                --input_format b \
                --iterations 30 \
                --weighting 3.0 \
                --header \
		--err_log $(DIR_P_ERRORS)/kitsch.$(*F).errors \
                < $*.sub_matrix \
		> $*.kitsch
	@touch $*.kitsch_finished
	@echo "	kitsch finish"

# Rejoin Kitsch results
$(DIR_P_TREE)/subbranches.rejoined: $(KITSCH_RESULTS_FILES)
	@echo "	rejoin_subbranches submitted"
	@for dir in $(DIR_P_TREE)/*.kitsch; do echo $$dir; done | xargs cat > $(DIR_P_TREE)/rejoined.all_kitsch
	@$(DIR_BIN)/tree_rejoin_subbranches \
		 --output_path "$(DIR_P_TREE)/[FILENAME].rejoined" \
		 --err_log $(DIR_P_ERRORS)/rejoin_subbranches.errors \
		 < $(DIR_P_TREE)/rejoined.all_kitsch
	@date > $(DIR_P_SUMMARY)/kaks.upgma_kitsch_finish.date
	@touch $(DIR_P_TREE)/subbranches.rejoined
	@#consolidate upgma error files
	$(call consolidate_non_empty_files,$(DIR_P_ERRORS),extract_sub_tree.clusters.*.errors,extract_sub_tree.errors)
	@#consolidate upgma error files
	$(call consolidate_non_empty_files,$(DIR_P_ERRORS),upgma.clusters.*.errors,upgma.errors)
	@#consolidate kitsch error files
	$(call consolidate_non_empty_files,$(DIR_P_ERRORS),kitsch.clusters.*.errors,kitsch.errors)
	@echo "	rejoin_subbranches fragments finish"
	@echo

# Interpret all rejoined / unchanged results in terms of orthology
$(DIR_OUTPUT_INTERPRET_TREE)/ds.tree: $(DIR_P_TREE)/subbranches.rejoined
	@echo "	copying trees to $(DIR_OUTPUT_INTERPRET_TREE)/ds.tree"
	@find $(DIR_P_TREE) -name *.rejoined | xargs cat > $@ 
	@find $(DIR_P_TREE) -name *.unchanged | xargs cat >> $@ 
	@echo "	Finished copying trees"

# Interpret all rejoined / unchanged results in terms of orthology
$(DIR_OUTPUT_INTERPRET_TREE)/transcript.ortholog.clades: $(DIR_OUTPUT_INTERPRET_TREE)/ds.tree
	@echo "	interpret trees"
	@echo 
	@$(DIR_BIN)/tree_interpret_nh \
		--x_taxon_regex "($(REGEX1)[^ \t]+)" \
		--y_taxon_regex "($(REGEX2)[^ \t]+)" \
		--output_path "$(DIR_OUTPUT_INTERPRET_TREE)/[FILENAME].results" \
		--x_taxon_name $(SPECIES1) \
		--y_taxon_name $(SPECIES2) \
		--orphaned_clades $(DIR_OUTPUT_INTERPRET_TREE)/orphaned.clades \
		--err_log $(DIR_P_ERRORS)/interpret.orthology.errors \
		--verbose \
		< $< \
		> $@
	@echo "	interpret trees finish"
	@echo


# Interpret all rejoined / unchanged results in terms of orthology
$(DIR_P_KAKS_RESULTS)/all.ortholog_ids: $(DIR_OUTPUT_INTERPRET_TREE)/transcript.ortholog.clades
	@echo "	get all identifier pairs for orthologs"
	@$(DIR_PROJECTBIN)/ortholog_sets_to_id_pairs.pl \
                        $< > $@
	@echo "	Done"
	@echo


$(DIR_OUTPUT)/ds_res.orthologs.both_filtered: $(DIR_P_KAKS_RESULTS)/all.ortholog_ids $(DIR_OUTPUT)/ds_res.both_filtered
	@echo "	get Ks for orthologs"
	@$(DIR_BIN)/seq_pairs_matching \
		--file1in $(DIR_OUTPUT)/ds_res.both_filtered \
		--file2in $< \
		--shared_in_full $@ \
		--file2_only_ids $(DIR_P_ERRORS)/orthologs.without.kaks_results.errors \
		--verbose
	@echo "	Done"
	@echo






#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#

#	Gene counts

#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#-----------------------------------------------------------------------------------------
#
#
# Filter by dS using 75% of shorter and 65% of longer and dS < 2.5
#
$(DIR_OUTPUT)/ds_res.2_5.both_filtered: $(DIR_OUTPUT)/dnds_res.both_filtered
	@echo "	Filter dS : 65% coverage of longer, 75% of shorter sequence ds < 2.5 ..."
	@$(DIR_BIN)/seq_pairs_filter_kaks_results \
                --ds_threshold 0.25 \
                --err_log $(DIR_P_ERRORS)/kaks.filter.both.errors \
		--coverage_residues 0.6 \
		--coverage_fragment 0.75 \
		--verbose \
		--output_ds_only \
		--results_summary $(DIR_P_SUMMARY)/kaks.filter.both.counts \
		--input_file $< \
		> $@
	@echo "	Filter Kaks finished"
	@echo


#_________________________________________________________________________________________
#
# All genes with blast results
#
$(DIR_P_GCNT)/1.blast_filtered.gene_counts.summary: $(DIR_P_BLASTRES)/curr.blastall_results $(DIR_OUTPUT)/ds_res.orthologs.both_filtered

	@echo "	Count all genes with blast results with significant overlap..."
	@$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_P_SEQ12)/cdna_hash_to_gene_id_species.map \
                --err_log $(DIR_P_ERRORS)/gene_count.blast_filtered.errors \
                --gene_output_path $(DIR_OUTPUT_GENE_LISTS)/blast_filtered.[SPECIES].genes \
		--verbose \
		--results_summary $@ \
		--input_file $<
	@echo



#_________________________________________________________________________________________
#
#  filtered kaks results
#
$(DIR_P_GCNT)/2.kaks_filtered.gene_counts.summary: $(DIR_OUTPUT)/ds_res.2_5.both_filtered $(DIR_P_GCNT)/1.blast_filtered.gene_counts.summary
	@echo "	Count all genes with filtered kaks results ..."
	@$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_OUTPUT)/prot_to_gene_ids.map \
		--err_log $(DIR_P_ERRORS)/gene_count.kaks_filtered.errors \
		--verbose \
		--gene_output_path $(DIR_OUTPUT_GENE_LISTS)/kaks_filtered.[SPECIES].genes \
		--results_summary $@ \
		--input_file $<
	@echo

#_________________________________________________________________________________________
#
#  orthologs
#
$(DIR_P_GCNT)/3.ortholog_clades.gene_counts.summary: $(DIR_P_KAKS_RESULTS)/all.ortholog_ids $(DIR_P_GCNT)/2.kaks_filtered.gene_counts.summary
	@echo "	Count all genes with ortholog clade results ..."
	@$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_OUTPUT)/prot_to_gene_ids.map \
		--gene_output_path $(DIR_OUTPUT_GENE_LISTS)/ortholog_clades.[SPECIES].genes \
		--err_log $(DIR_P_ERRORS)/gene_count.ortholog_clades.errors \
		--verbose \
		--results_summary $@ \
		--input_file $<
	@echo



#_________________________________________________________________________________________
#
# All genes with blast results in the other species
#
$(DIR_P_GCNT)/4.blast_other_species.gene_counts.summary: $(DIR_P_BLASTRES)/curr.blastall_results $(DIR_P_GCNT)/3.ortholog_clades.gene_counts.summary
	@echo "	Count all genes with significant blast results in the other species..."
	@$(DIR_BIN)/seq_pairs_substitute_identifiers \
		--input_file $< \
		--map_identifiers $(DIR_P_SEQ12)/cdna_hash_to_prot_id.map \
                --err_log $(DIR_P_ERRORS)/gene_count.blast_other_species_map.errors | \
	$(DIR_BIN)/seq_pairs_filter_by_regex \
                --err_log $(DIR_P_ERRORS)/gene_count.blast_other_species_ortho_filter.errors \
		--regex "($(REGEX1)|$(REGEX2))" \
		--patterns "$(REGEX1)" \
		--alt_patterns "$(REGEX2)"  | \
	$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_OUTPUT)/prot_to_gene_ids.map \
                --err_log $(DIR_P_ERRORS)/gene_count.blast_other_species_list.errors \
                --gene_output_path $(DIR_OUTPUT_GENE_LISTS)/blast_other_species.[SPECIES].genes \
		--verbose \
		--results_summary $@
	@echo




#_________________________________________________________________________________________
#
#  filtered kaks results in the other species
#
$(DIR_P_GCNT)/5.kaks_filtered_other_species.gene_counts.summary: $(DIR_OUTPUT)/ds_res.2_5.both_filtered $(DIR_P_GCNT)/4.blast_other_species.gene_counts.summary
	@echo "	Count all genes with filtered kaks results in the other species..."
	@$(DIR_BIN)/seq_pairs_filter_by_regex \
		--input_file $< \
                --err_log $(DIR_P_ERRORS)/gene_count.kaks_filtered_other_species_ortho_filter.errors \
		--regex "($(REGEX1)|$(REGEX2))" \
		--patterns "$(REGEX1)" \
		--alt_patterns "$(REGEX2)"  | \
	$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_OUTPUT)/prot_to_gene_ids.map \
		--err_log $(DIR_P_ERRORS)/gene_count.kaks_filtered_other_species.errors \
		--verbose \
		--gene_output_path $(DIR_OUTPUT_GENE_LISTS)/kaks_filtered_other_species.[SPECIES].genes \
		--results_summary $@
	@echo

#_________________________________________________________________________________________
#
#  orthologs
#
$(DIR_P_GCNT)/6.ortholog_clades_other_species.gene_counts.summary: $(DIR_P_KAKS_RESULTS)/all.ortholog_ids $(DIR_P_GCNT)/5.kaks_filtered_other_species.gene_counts.summary
	@echo "	Count all genes with orthologs results in the other species ..."
	@$(DIR_BIN)/seq_pairs_filter_by_regex \
		--input_file $< \
                --err_log $(DIR_P_ERRORS)/gene_count.ortholog_clades_other_species_ortho_filter.errors \
		--regex "($(REGEX1)|$(REGEX2))" \
		--patterns "$(REGEX1)" \
		--alt_patterns "$(REGEX2)"  | \
	$(DIR_BIN)/seq_pairs_list_identifiers \
		--map_prot_to_gene_to_species $(DIR_OUTPUT)/prot_to_gene_ids.map \
		--gene_output_path $(DIR_OUTPUT_GENE_LISTS)/ortholog_clades_other_species.[SPECIES].genes \
		--err_log $(DIR_P_ERRORS)/gene_count.ortholog_clades_other_species.errors \
		--verbose \
		--results_summary $@
	@echo


#_________________________________________________________________________________________
#
#  Counts summary
#
$(DIR_OUTPUT)/phyop_stages.gene_counts.summary:  $(DIR_P_GCNT)/6.ortholog_clades_other_species.gene_counts.summary
	@echo "	Summarise all counts to $@ ..."
	@$(CURDIR)/parse_gene_counts.pl $(DIR_P_GCNT)/*.gene_counts.summary > $@
	@cp $@ $(DIR_P_SUMMARY)/phyop_stages.gene_counts.summary
	@echo "	Done"
	@echo
	


#
#   OUTPUT files
#
fifth_part: $(DIR_OUTPUT)/phyop_stages.gene_counts.summary \
            $(DIR_OUTPUT)/ds_res.2_5.both_filtered \
	    $(DIR_OUTPUT)/dnds_res.both_filtered \
	    $(DIR_OUTPUT_INTERPRET_TREE)/ds.tree \
            $(DIR_OUTPUT_INTERPRET_TREE)/transcript.ortholog.clades \
            $(DIR_OUTPUT)/ds_res.orthologs.both_filtered
#
# End of fifth part
#
#-----------------------------------------------------------------------------------------

